1 Area of Interest

2 Landsat Tiles

The dimensions of the stacked image are 7,811 rows, 7,861 columns, 59,996,291 cells and 6 layers. The CRS is UTM on the WGS84 datum. The resolution is 30m x 30m.

The dimensions of the stacked image are 340 rows, 346 columns, 117,640 cells and 6 layers. The CRS is UTM on the WGS84 datum. The resolution is 30m x 30m.

3 Plots

3.1 Non-Contrasted

3.1.1 Natural Color

3.1.2 Color Infared

3.1.3 False Color Water Focus

3.1.4 Land From Water Focus

3.2 Contrasted

3.2.1 Natural Color

3.2.2 Color Infared

3.2.3 False Color Water Focus

3.2.4 Land From Water Focus

4 Flooding

4.1 Raster Algebra

A we can see, the NWDI, MNDWI, and WRI are all very similar in the way they pick up water cells with the same color, but they deviate differently in the dry ground cells. NDVI picks up water in a different color, and land in the opposite as well. The SWI seems to pick up only water cells and no land cells.

4.2 Thresholding

SWI has multiple NA values so we will have to replace those.

5 Clustering

This shoes us they are extracted as a matrix of size 117640 x 6

##    
##         1     2     3     4     5     6
##   0   488 17474  8361 53172 22289   655
##   1 13387     0    50   108  1646    10

6 Flooding Stats

Flood Area for Different Methods
Flood Area m^2
NDVI 5999400
NDWI 6490800
MNDWI 10745100
WRI 7622100
SWI 13680900
K.Means 12487500

6.1 Flood Uncertainty

6.1.1 Base R

6.1.2 Open Steet Map

Some of the values in the raster are not whole numbers because the color palette is continuous, so to make the overall raster, which is discrete, fit this color scheme, it is possible it combines values nearby and makes a raster value a non-whole number, which results in the color smoothing of the map.

7 Appendix

knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)
library(tidyverse)
library(sf)
library(raster)
library(getlandsat)
library(mapview)
library(osmdata)
library(knitr)
bb = read_csv('../data/uscities.csv') %>%
  filter(city == 'Palo') %>%
  st_as_sf(coords = c('lng','lat'), crs = 4326) %>%
  st_transform(5070) %>%
  st_buffer(5000) %>%
  st_bbox() %>%
  st_as_sfc()
meta = read_csv('../data/palo-flood.csv')

files = lsat_scene_files(meta$download_url) %>%
  filter(grepl(paste0('B',1:6,'.TIF$', collapse = '|'), file)) %>%
  arrange(file) %>% 
  pull(file)

st = sapply(files, lsat_image) 

b = stack(st) %>% setNames(paste0('band', 1:6))
cropper = bb %>% st_as_sf() %>% st_transform(crs(b))

cr_b = crop(b, cropper)

r <- cr_b %>% setNames(c('coastal', 'blue', 'green', 'red', 'nir', 'swir'))

plotRGB(r, r = 4, g= 3, b = 2)
plotRGB(r, r = 5, g= 4, b = 2)
plotRGB(r, r = 5, g= 6, b = 4)
plotRGB(r, r = 7, g= 5, b = 3)
plotRGB(r, r = 4, g= 3, b = 2, stretch = 'hist')
plotRGB(r, r = 5, g= 4, b = 2, stretch = 'hist')
plotRGB(r, r = 5, g= 6, b = 4, stretch = 'lin')
plotRGB(r, r = 7, g= 5, b = 3, stretch = 'lin')
ndvi = (r$nir - r$red)/(r$nir + r$red)

ndwi = (r$green - r$nir)/(r$green + r$nir)

mndwi = (r$green - r$swir)/(r$green + r$swir)

wri = (r$green + r$red)/(r$nir + r$swir)

swi = 1/(sqrt(r$blue-r$swir))

rs <- stack(ndvi, ndwi, mndwi, wri, swi) %>% setNames(c('NDVI', 'NDWI', 'MNDWI', 'WRI', 'SWI'))

plot(rs, col = colorRampPalette(c("blue", "white", "red"))(256))

ndvi_t = function(x) {ifelse(x<=0,1,0)}
ndwi_t = function(x) {ifelse(x>=0,1,0)}
mndwi_t = function(x) {ifelse(x>=0,1,0)}
wri_t = function(x) {ifelse(x>=1,1,0)}
swi_t = function(x) {ifelse(x<=5,1,0)}


ndvi_f = calc(ndvi, ndvi_t)
ndwi_f = calc(ndwi, ndwi_t)
mndwi_f = calc(mndwi, mndwi_t)
wri_f = calc(wri, wri_t)
swi_f = calc(swi, swi_t)

f_stack = stack(c(ndvi_f, ndwi_f, mndwi_f, wri_f, swi_f)) %>% setNames(c('NDVI', 'NDWI', 'MNDWI', 'WRI', 'SWI'))
plot(f_stack, col = colorRampPalette(c("white","blue"))(256))
#sum(is.na(values(ndvi_f)))
#sum(is.na(values(ndwi_f)))
#sum(is.na(values(mndwi_f)))
#sum(is.na(values(wri_f)))
#sum(is.na(values(swi_f)))
values(swi_f) <- ifelse(is.na(values(swi_f)), 0, values(swi_f) )

#sum(is.na(values(swi_f)))
vals <- getValues(r)

#dim(vals)
vals <- na.omit(vals)
vs <- scale(vals)
idx <- 1:dim(vals)[1]
set.seed(20200905)

rast.clust <- kmeans(vs, 6, iter.max = 100)

kmeans_raster <- rs$MNDWI
values(kmeans_raster) <- rast.clust$cluster

plot(kmeans_raster)
table(getValues(swi_f), getValues(kmeans_raster))
k_thresh = function(x) {ifelse(x==1,1,0)}

k_flood = calc(kmeans_raster, k_thresh) %>% setNames('K-Means')

f_stack <- addLayer(f_stack, k_flood)
plot(f_stack, col = colorRampPalette(c("white","blue"))(256))
flood_area <- cellStats(f_stack, 'sum') * (30^2)
flood_area_km <- flood_area*0.000001

#names(flood_area)

kable(data.frame(flood_area), col.names = 'Flood Area m^2', caption = 'Flood Area for Different Methods')
rast_uncert <- calc(f_stack, function(x){sum(x)}) %>% setNames('Flood Uncertainty')
plot(rast_uncert, col = blues9, main = 'Flood Uncertainty')
values(rast_uncert) <- ifelse(values(rast_uncert)==0, NA, values(rast_uncert))
mapview(rast_uncert)